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It has been recently shown numerically that the transition from integrability to chaos in quantum 
systems and the corresponding spectral fluctuations are characterized by noise with I < a < 2. 
The system of interacting trapped bosons is inhomogeneous and a complex system. The presence 
of external harmonic trap makes it more interesting as in the atomic trap the bosons occupy partly 
degenerate single-particle states. Earlier theoretical and experimental results show that at zero 
temperature the low-lying levels are of collective nature and high-lying excitations are of single 
particle nature. We observe that for few bosons, P(s) distribution shows the Shnirelman peak which 
exhibits a large number of quasi-degenerate states. For large number of bosons the low-lying levels 
are strongly affected by the interatomic interaction and the corresponding level fluctuation shows a 
transition to Wigner with increase in particle number. It does not follow GOE (Gaussian Orthogonal 
Ensemble) Random Matrix predictions. For high-lying levels we observe the uncorrelated Poisson 
distribution. Thus it may be a very realistic system to prove that i noise is ubiquitous in nature. 

PACS numbers: 03.75.Hh, 31.15.Ja, 05.45.Mt, 05.45.Pq, 05.45.Tp 



I. INTRODUCTION 

Although there is no precise definition of quantum 
chaos, however it is closely related with the energy level 
fluctuation properties of a quantum system. Bohigas 
conjectured that level fluctuation of quantum system 
whose classical limit is chaotic, is described by the 
random matrix theory (RMT) whereas spectral 

fluctuation of classically intcgrable system obeys Poisson 
statistics @- The concept of quantum chaos plays an 
important role in the understanding of the universal 
properties of the energy level spectrum of quantum 
system. However the complex natural systems are 
neither fully intcgrable, nor fully chaotic and attains 
special interest. RMT introduced by Wigner has been 
widely used in the description of complex spectrum of 
atomic nucleus, atoms and molecules [J-Q. On the other 
hand, the bosonic ensembles in the dense limit may be 
ergodic with increase in the number of single particle 
states @. In last few years, interacting bosonic systems 
got special interest due to the experimental observation 
of Bose-Einstein condensation |7HlCf. The presence 
of external harmonic trap makes it more interesting as 
stated by Asaga that in an atomic trap, bosonic atoms 
occupy partly degenerate single-particle states [ll|. Al- 
though it is argued that random-matrix approach should 
reveal the generic features of the spectrum however there 
is neither analytical treatment, nor systematic numerical 
calculations in this direction. The chaotic signature 
in the time evolution of Bosc-Einstcin condensation 
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driven by the time-periodic harmonic or kicked pulses 
are observed [12Hl4l |. But energy level statistics of the 
experimentally dilute BEC has not been studied yet. In 
the earlier analysis of Bohigas in nuclear and atomic 
spectra the nearest-neighbor spacing distribution agrees 
very well with GOE (Gaussian Orthogonal Ensem- 
ble) [l5j]. However for the interacting trapped bosons, 
it seems to contradict the usual expectation based on 
RMT. Very recently spectral properties of trapped ID 
ultracold fermions in optical lattices are studied and 
the interplay of repulsive interaction with the external 
harmonic trap is observed [l6|. So it is also very much 
interesting to study the energy-level statistics of trapped 
bosons which are spatially inhomogeneous and we may 
expect new and rich physics. 

Recently a different approach to characterize quan- 
tum chaos has been proposed based on the idea that the 
corresponding energy level sequence is analogous to the 
discrete time series. The level fluctuation is well charac- 
terized by the Fourier power spectrum and a power law 
behavior has been identified. This is conjectured that 
spectral fluctuations of chaotic quantum systems are 
characterized by j noise whereas complete integrable 
systems exhibit -p noise fl7l - [20l ]. The earlier studies in 
this direction involve quantum billiards, nonintegrable 
coupled quartic oscillator, kicked top, integrable spin 
chain 0, I2ll - l23l |. In this present work we study the 
system of N interacting bosons at zero temperature in 
the presence of external trap. The choice of such system 
is important for various reasons. Firstly: it is inhomo- 
geneous and complex system due to the presence of two 
energy scales. Interatomic interaction is characterized 
by Na s , where a s is the s-wave scattering length and the 
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external trap energy is characterized by fiw, where lo is 
the external trap frequency. From the earlier theoretical 
and experimental results it is an established fact that at 
zero temperature the low-lying collective excitations are 
strongly affected by the interatomic interaction when 
the high-lying excitations are of single particle nature 
p2l [2J-[2g. The transition from collective to single 
particle excitations makes us more curious to study 
the level fluctuation and to verify whether noise is 
ubiquitous in nature. Secondly: the system directly 
manisfests the experimental Bosc-Einsticn condensation 
0-tHj- For the present calculation we consider the 
iV-body bosonic system at zero temperature. There 
may be a very small effect of thermal cloud around the 
condensate even at zero temperature and the condensate 
is depleted due to the interaction [27} • However for 
the present calculation we ignore that as the whole 
condensate is described by a single and fixed scattering 
length and the condensate is extremely dilute. Thus 
the effect of damping does not appear in our present 
calculation. However the effect of damping may be 
important when the interaction is tuned by external 
magnetic field. Thus the system in our present work 
is neither fully chaotic (for low-lying levels) nor fully 
integrable (for high-lying levels) due to the interplay of 
two-energy scales. At this point we should mention that 
Bohigas analysed thoroughly the nuclear shell-model 
and neutron resonance data for different nuclei. The 
nearest-neighbor spacing distribution of the nuclear 
data ensemble (NDE) agrees very well with the GOE 
prediction [l5[ . In the atomic spectra the levels with the 
same quantum numbers also show Wigner type spacing 
distribution. Thus in nuclear and atomic spectra, the 
regular features of low-lying part of the spectrum and 
chaotic features of high-lying collective levels are well 
established fact. However for the interacting trapped 
bosons, it seems to contradict the usual expectation 
based on RMT as for the experimental BEC, the 
low-lying excitations are collective where the interatomic 
interaction plays a crucial role and the high-lying levels 
are of single particle nature due to the dominating effect 
of external harmonic trap. 



The paper is organized as follows. Sec. II deals 
with the methodology which contains the many-body 
technique to calculate the energy levels. Choice of 
interaction and the correlation function are also dis- 
cussed in the same section. In the Sec. Ill we discuss 
several statistical tools and results. Sec. IV concludes the 
summary. 



II. METHODOLOGY 

A. Many-body calculation with potential harmonic 

basis 

In order to calculate the energy levels of the con- 
densate we solve the Schrodinger equation by our 
newly developed correlated potential harmonic expan- 
sion method (CPHEM) with a short-range correlation 
function. CPHEM has already been established as a 
very successful technique for the study of dilute BEC 
[28M301 ]. In this method we keep all possible two-body 
correlation and also use a realistic interatomic interac- 
tion which is clearly an improvement over the mean-field 
Gross-Pitaevskii (GP) theory [U,|3l|. We briefly discuss 
the technique below. 

We consider a system of A = (N + 1) identical bosons 
interacting via two-body potential V{fij) = V(?i — fj) 
and confined in an external harmonic potential of fre- 
quency oj. The time- independent quantum many-body 
Schrodinger equation is given by 

-E\*(r u ...,f A )=0; (1) 

where m is the mass of the each boson and E is the 
energy of the condensate. After eliminating the center of 
mass motion by using the standard Jacobi vectors [32I — 
H3, defined by 



cwdrr^-iE^) (< = l •••*). ( 2 ) 

we obtain the relative motion of TV-body system as 




1=1 



-E R ]*(( U ...,( N ) = 0, (3) 

Vt ra p is the effective external trapping potential, Vint is 
the sum of all pair-wise interactions expressed in terms 
of the Jacobi vectors and En is the relative energy of the 
system i.e. E = En + §7iw. 

Now it is to be noted that Hyperspherical harmonic ex- 
pansion method (HHEM) is an ab-initio tool to solve the 
many-body Schrodinger equation where the total wave 
function is expanded in the complete set of hyperspheri- 
cal basis [32[ . Although HHEM is a complete many-body 
approach and includes all possible correlations, it can not 
be applied to a typical BEC containing few thousands to 
few millions of bosons. Due to the large degeneracy of the 
HH basis, HHEM is restricted only to three-particle sys- 
tems [32|, |35( . Since the typical experimental BEC is de- 
signed to be very dilute and the probability of three and 
higher-body collisions is negligible, we can safely ignore 
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the effect of three and higher-body correlations. There- 
fore only two-body correlation and pairwise interaction 
among the bosons is important. It allows us to decom- 
pose the total wave function ^ into two-body Faddeev 
component for the interacting (ij) pair as 



* = E <M% r ) 



(4) 



It is worth to note that 0y is a function of two-body 
separation (fi,j ) and the global hyperradius r is given by, 

r = Q- Thus the effect of two-body correlation 

comes through the two-body interaction in the expansion 
basis. <f>ij is symmetric under P^- for bosons and satisfy 
the Faddeev equation 



[T + Vtrap - E R ] <P l3 = -V(f 



A 

E ■ 

k,l>k 



(5) 



A' 



where T = — ^ J^V 2 -. is the total kinetic energy. Oper- 

i=l 

ating J2i j>i on both sides of equation (5), we get back 
the original Schrodingcr equation. In this approach, we 
assume that when (ij) pair interacts, the rest of the 
bosons are inert spectators. Thus the total hyperangular 
momentum quantum number as also the orbital angu- 
lar momentum of the whole system is contributed by the 
interacting pair only. Next we expand <f>ij in the sub- 
set of hypcrspherical harmonics (HH) necessary for the 
expansion of V{fij). 



non-interacting pair 



h — h — ■■■ — In-i — 0, 
mi = ni2 = ... = toa-i = 0, 
n 2 = n 3 = ...njv-i = 0, 



(9) 
(10) 
(11) 



and for the interacting pair ijv = I, tojv = m and 
tin = K. Thus the 3N dimensional Schrbdinger equation 
reduces effectively to a four dimensional equation with 
the relevant set of quantum numbers: hyperradius r, or- 
bital angular momentum quantum number I, azimuthal 
quantum number m and grand orbital quantum number 
2K+1 for any N. Substituting Eq(6) into Eq.(5) and pro- 
jecting on a particular PH, a set of coupled differential 
equation (CDE) for the partial wave u l K (r) is obtained 



m dr 2 



V trap (r) + ^{C(C + 1) 
+4K(K + a + /3 + 1)} - E R ] U K i(r) 



(12) 



K 



fKiV KK >(r)fK>iU K >i(r) = 



U K i = fi<iu l K {r), a 



3A-8 



and 



where C = I + ^ 
= 1 + 1/2. 

Jki is a constant and represents the overlap of the PH 
for interacting partition with the sum of PHs correspond- 
ing to all partitions (34|. The potential matrix element 
Vkk> (r) is given by 



2K' + 1 



N )dQ N 



(13) 



(6) 



K 



denotes 



the full set of hyperangles in the 3N- 
dimensional space corresponding to the (ij)th interacting 
pair and ^ir+lC^Jv) * s cadled the PH basis. It has an an- 
alytic expression: 



V. 



l,m 
2K+l\ 



l 2 i +l (ct>)y (D-3); 



D = 3N, 
(7) 

yo(D — 3) is the HH of order zero in the (3N — 3) dimen- 
sional space spanned by {d, ...,Ca-i} Jacobi vectors; <j> 
is the hyperangle given by = r coscj). For the remain- 
ing (N — 1) noninteracting bosons we define hyperradius 



Pij 



\ 



N-l 

E 



(-2 



(8) 



= rfj +Pij and r represents the global hyper- 



such that r 2 

radius of the condensate. The set of (3N — 1) quantum 
numbers of HH is now reduced to only 3 as for the [N— 1) 



B. Choice of interaction and introduction of 
additional short range correlation 

In the mean-field GP equation the two-body inter- 
action is taken as the contact 5 potential, the interac- 
tion strength being proportional to the s-wave scattering 
length a s . A positive value of a s gives a repulsive conden- 
sate and a negative value of a s gives an attractive conden- 
sate. But the contact interaction completely disregards 
the detailed structure. However a realistic interatomic 
interaction, like the van der Waal potential, is always as- 
sociated with an attractive — tail at large separation 
and a strong repulsion at short separation. Depending on 
the nature of these two parts, a s can be either positive or 
negative. In our earlier calculations (36j we have already 
observed the effect of shape-dependent interatomic inter- 
action in the many-body calculation. So for our present 
calculation we choose the van der Waal potential with a 
hard core repulsion of radius r c , viz, Vfaj) = oo for r. 



< r c and —j^p for nj > r c . The value of Cq is fixed for 

a given system and for 87 Rb atoms Cq — 6.4898 x 10~ n 
o.u. [31[. Throughout our calculation we choose aho 
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as the unit of length (o.u.) and energy is also 

expressed in the unit of oscillator energy (hoj). For a 
given two-body interaction a s can be obtained from the 
solution of two-body equation with zero energy. 



H 2 1 d 



r l^l ]+V{rijHnj) = . (14 ) 



13 dr t 



The solution of the two-body equation shows that the 
value of a s changes from negative to positive and thus 
passing through an infinite discontinuity as r c decreases 
[Fig. 1]. At each discontinuity one extra node appears 
in the two-body wave function which corresponds to one 
extra two-body bound state. With a tiny increase in 
r c , across the infinite discontinuity a s changes drastically 
from a very large positive value to a large negative value 
and the properties of the condensate changes drastically 
[3~fl ] . In the GP equation one uses a s directly without any 
such detailed knowledge of actual interatomic potential. 
For the present calculation we choose a s =0. 00433 o.u. 
which mimics the JILA experiment with 87 Rb atoms Q . 
The corresponding value of r c is 1.121 xl0~ 3 o.u. which 
causes one node in the two-body wave function. The nor- 
malization constant is chosen to make the wave function 
positive at large r^. 

In the experimental BEC, the Bose gas is extremely 
dilute, the average interparticle separation is much larger 
than the range of the two-body interaction. This is re- 
quired to prevent the three-body collision and formation 
of molecules. Thus the pair of particles with practically 
zero kinetic energy do not come closer than a s . Whereas 
the zero th order PH is a constant [33[ and will give a 
large probability even for ry— > 0, it causes very slow con- 
vergence in the PH basis [Eq.(6)]. To compensate this 
we additionally include a short range correlation func- 
tion r](rij) in the PH expansion. As the fundamental 
assumption in our method is to consider only (ij) pair 
interaction when the remaining particles are simply in- 
ert spectators, the correlation function is obtained as the 
zero-energy solution of the two-body equation [Eq.(14)]. 
The correlation function quickly attains asymptotic form 
(1 — ^r) for large r^. We replace Eq(6) by 



J2 V ^ + i(^W K (r) V ( ril ) ■ (15) 



K 



The correlated PH (CPH) basis becomes 



correlated 



2K 



^H^), (16) 



The correlated potential matrix Vkk> ( r ) is now given by 
V KKI (r) = (hfh a f f,)-ix 

j^{Pf(z)V (r^) Pg(z)r, (ryf^z) W,(z)}dz ■ 

(17) 

Here P^(z) is the Jacobi polynomial, and its norm and 
weight function are and W\{z) respectively (37| . 



One may note that the inclusion of f](rij) makes 
the PH basis non-orthogonal. One may surely use the 
standard procedure for handling non-orthogonal basis. 
However in the present calculation we have checked that 
r)(rij) differs from a constant value only by small amount 



and the overlap 

quite small. Thus we get back the Eq(12) approximately 
when the correlated potential matrix is calculated by 
Eq(17). 

Finally the coupled differential equation (CDE), 
Eq. (12), is solved by the hyperspherical adiabatic ap- 
proximation (HAA) [23| . In HAA, one assumes that the 
hyperradial motion is slow compared to the hyperangu- 
lar motion. Hence the latter is separated adiabatically 
and solved for a particular value of r, by diagonalizing 
the potential matrix together with the diagonal hyper- 
centrifugal repulsion in Eq. (12). The lowest eigenvalue, 
a>o(r) is the effective potential for the hyperradial mo- 
tion and in this effective potential the entire condensate 
moves as a single entity. Thus in HAA, the approximate 
solution (the energy and wave function) of the conden- 
sate is obtained by solving a single uncoupled differential 
equation 



h2 d2 / x „ 

TT + w oW - E R 

m dr z 



<o(r) = , 



(18) 



subject to appropriate boundary conditions on Co(^)- 
The function Co (r) is the collective wave function of the 
condensate in the hyperradial space. The lowest lying 
state in the effective potential u>o(r) corresponds to the 
ground state of the condensate. The total energy of 
the condensate is obtained by adding the energy of the 
center of mass motion (§fiw) to Er. 

Thus by employing the CPHEM and HAA we re- 
duce the multi-dimensional problem into an effective 
one-dimensional problem in hyperradial space and the 
effective potential u>o(r) provides both the qualitative 
and quantitative description of the system. As in 
our many-body picture, the collective motion of the 
condensate is characterised by the effective potential, 
the excited states in this potential are the states with 
n th radial excitation and I th surface mode and arc 
generally denoted by E n \. Thus -Eoo corresponds to the 
ground state and / ^ corresponds to several surface 
modes. For / > 0, we calculate the potential matrix from 
the diagonal hypercentrifugal term. We have checked 
that the contribution coming from the off-diagonal 
matrix element is very very small and we disregard 
these matrix elements as they make the computation 
very slow. The calculation of low-lying collective modes 
are in good agreement with the experimental results 
and other calculations [H, Ho|. For energy much 
larger than the chemical potential (fi) we observe that 
the states are separated at energy close to harmonic 
oscillator energies (~ Hu>). This transition from the 
low-energy collective modes to high-lying single particle 
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excitation arc further used for the statistical calculations. 



III. RESULTS 

The integrated level density N(E) has two parts. 
One is the smooth part (N(E)) and a fluctuating part 
(N(E)). To compare the fluctuation of different systems 
or different parts of the same system, the smooth part is 
removed by the unfolding procedure. Unfolding maps the 
energy levels Ei to e.; with the unit mean level density. 
For the present analysis the many-body level density is 
approximated by a polynomial and unfolding is done by 
7 th order polynomial. Wc unfold each spectrum sepa- 
rately for a specific value of / and form an ensemble hav- 
ing the same symmetry. Then the nearest neighbor spac- 
ing is calculated as Sj = e.;+i — e,, i = 1, 2, ...n. For the 
further study of correlation and level-repulsion between 
energy levels we utilize the established analogy between 
the energy spectrum and discrete time series [17i - ll9l l4l| . 
The energy spectrum is considered as a discrete signal 
and the fluctuations of the excitation energy as discrete 
time series. The S n statistics has been used in RMT to 
study how the consecutive level spacings are correlated. 
It is defined as 

n 

S„ = 2j(sj- < s >) = e n+ i - ei - n (19) 

i=l 

As the average value of Si is < s >= 1, S n represents the 
deviation of (n + l) th level from the mean value i.e. the 
fluctuation of (n + l) th excited state. It is also closely 
related to the level density fluctuations and one can write 
S n =— N(E n ±i) if the ground state energy is shifted ap- 
propriately |4l| . Thus it represents the accumulated level 
density fluctuation at E=E n +x. S n is similar to the time 
series and n represents the discrete time [l8l - [20l l4lj . The 
power spectrum is then defined as the square modulus of 
the Fourier transform as 

^H^EW-— )l 2 (20) 

where k = 1,2, ....n and / = represents the fre- 
quency and M is the size of the series Therefore, 
the statistical behavior of level fluctuation can be 
established by < P% > statistic which measures both 
short and long range correlation. It is verified that 
the power laws Pf_oc4 both for fully chaotic and 
integrable systems |r7H20| . But depending on the level 
correlation in the chosen system a scales smoothly from 
1 (chaoti c sy stem) to 2 (for uncorrelated and integrable 
system) [17H20l | . However in the integrable spin chains 
of Halden-Shastry type, the spectral fluctuations exhibit 
■jj noise rather than the expected -p- noise [23[ . 

In Fig. 2 we display the energy level fluctuations 
for different number of energy levels for 5000 bosons in 



the trap. For the low-lying levels we expect level corre- 
lation. As the low-lying levels are highly affected by the 
interatomic interaction, the energy spectra shows level 
repulsion and strong spectral rigidity. This is reflected in 
the Fig. 2(a) which looks like the antipersistent time se- 
ries for the lowest 500 levels. The 5 n statistics for the low 
levels is very close to the GOE spectra which indicates 
high level correlation due to the interatomic interaction. 
For the intermediate levels, the effect of interatomic inter- 
action gradually decreases and the external trap starts to 
dominate. Thus the system is expected to show a mixed 
and complex statistics. When a part of the levels are 
correlated due to interatomic interaction and two-body 
correlation, the other part do not repel each other and 
uncorrelated. It is similar to the classical mixed system, 
where a part of phase space is completely regular with 
the other part chaotic. Thus the Fig. 2(b) shows that S n 
is neither persistent nor antipersistent. For much higher 
levels [Fig. 2(c)], the energy levels are uncorrelated due 
to the dominating effect of the external harmonic trap. 
The system is close to integrable and 8 n looks like a per- 
sistent series of Poisson spectra. To characterise long- 
range correlation in Fig. 2(d)-2(f), we plot the average 
values of the power spectrum < Pjt > for the same num- 
ber of levels as reported in Fig. 2(a)-2(c). It shows that 
the power spectrum follows the scaling law < P^ > ~ 
j^. The value of a is presented in Fig. 2(d)-2(f) for 
different number of levels. For low-lying correlated lev- 
els a = 1.31, for intermediate levels a = 1.72 and for 
high- lying levels a = 1.99. Thus a not only measures 
the chaoticity of the system but it measures the degree 
of integrability for complex systems. At this point we 
should mention that in a nice attempt the momentum 
distribution and temporal power spectra of nonzero tem- 
perature Bose-Einstein condensate are calculated using 
the Gross-Pitaevskii equation (42[. The temporal power 
spectra also shows i form where a = 2 — y (D is the 
dimension of space) [H, EH . Next to compare the re- 
sult with the most popular and well known statistics, we 
calculate the nearest neighbor spacing distribution P(s) 
and plot in Fig. 3 and in Fig. 6. In an earlier attempt 
in this direction we have reported some preliminary re- 
sults on level-spacing distribution P(s) J44j. We have 
shown that due to interatomic correlation the lower lev- 
els are strongly affected by the interaction, however the 
higher levels are uncorrelated. But our earlier results do 
not prove the Asaga's statement which says that in an 
atomic trap, bosonic atoms are in partly degenerate sin- 
gle particle states. This makes us very curious to study 
in details how the small interatomic interaction will act 
as a perturbation and lift the degeneracy. This needs 
further numerical analysis for varying number of bosons 
and with increase in number of levels. In Fig. 3 wc plot 
the P(s) distribution for the lowest 100 levels for differ- 
ent number of bosons. For TV = 3 with a s = 2.09 x 10~ 4 
o.u., the effective interaction Na s is 8.7 x 10~ 4 o.u. The 
system is very close to integrable as the effect of such 
small interaction is masked due to the effect of external 
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harmonic trap. At zero temperature the interaction en- 
ergy for N — 3 is almost negligible compared to the trap 
energy. Thus the small interaction acts as a very small 
perturbation and the exact degeneracy in the external 
3D harmonic trap is lifted and it results to the existence 
of large quasi-degenerate states. P(s) distribution ex- 
hibits ^-type peak called as Shnirelman peak. In the year 
1993, Shnirelman showed that for systems with time re- 
versal symmetry should exhibit such a (5-function peak 
near s = in the P(s) distribution. It is known as the 
Shnirelman peak. This peak appears due to the presence 
of symmetry and separating levels by symmetry, one will 
get back Poisson distribution. This indicates the pres- 
ence of bulk quasi-degenerate states in the level spacing 
distribution. In the first verification of Shnirelman theo- 
rem, Chirikov and Shcpelyansky studied the kicked rota- 
tor on a torus with time- reversal symmetry (45j . Later 
the theorem is verified in a more real physical quantum 
system. The Calogero-like three-body problem was stud- 
ied where the hidden continuous symmetry was broken 
by adding a three-body interaction term [461 ]. With fur- 
ther increase in N gradually in the trap, the lower levels 
show level-repulsion and the system smoothly changes 
to close to integrability to nonintcgrability. The corre- 
sponding P(s) distribution smoothly changes to Wigner 
like distribution with increase in N. Due to strong inter- 
atomic interaction the system becomes more correlated 
and show level-repulsion. 

In our present problem of trapped, interacting 
bosons, the exact degeneracy comes from the external 
harmonic trap. However due to the weak interatomic in- 
teraction, the effect of exact degeneracy is gradually lifted 
and it results to the quasi-degeneracy when the number 
of bosons in the trap is quite small. For better resolution 
of the Shnirelman peak appeared in Fig. 3 (with N=3) 
we plotted the same in Fig. 4 in finer details. A huge peak 
in the first bin of the histograms clearly demonstrates the 
existence of global quasi-degeneracy in accordance with 
the Shnirelman theorem. In the top left-most panel in 
Fig. 4, we observe the peak has a finite width which is 
further associated with the Poissonian tail. This peak 
contains important information about the structure of 
the quantum system. The resolution of the peak is fur- 
ther plotted in Fig. 5 where we present the integral level- 
spacing distribution I(s) = NP(s), normalized to unity. 
It has two separate regions. The rightmost steep-increase 
of I(s) corresponds to the Poissonian tail of Fig. 4. The 
leftmost part is more interesting. It shows the linear 
dependence between / and ln(s), which represents the 
structure of the Shnirelman peak. 

The results for higher levels close to 4000 levels and 
for the same set of N values as reported in Fig. 3, are 
plotted in Fig. 6. For N — 3 , P(s) distribution again 
shows the sharp peak as expected. As the high lying ex- 
citations are of single particle nature, the energy levels 
are now uncorrelated and the corresponding P(s) dis- 
tribution shows Poisson type fluctuation with increase 
in number of bosons. It confirms that for higher lev- 



els the system again becomes close to intcgrable as the 
effect of external trap strongly dominates. The observa- 
tion is in correlation with the earlier observation of 8 n 
statistics and power spectrum. P(s) measures the short- 
range correlation. The A3 statistic is usually used to 
investigate the long-range correlation. It gives the sta- 
tistical measure of the rigidity of finite spectral level se- 
quence. For a given energy interval L, it is determined by 
the least square deviation of the staircase from the best 
straight line fits it. In Fig. 7 we plot the spectral average 
< A^(L) > for different energy levels. For higher energy 
levels < A^(L) > bends to Poisson whereas for low- lying 
collective levels it close to GOE prediction. 



IV. CONCLUSIONS 

In summary, we have shown that the analogy 
between quantum energy spectra and time series is an 
efficient and powerful way to characterize quantum level 
fluctuation. Although the statistical behaviour of level 
fluctuation and the corresponding power spectrum are 
understood for fully chaotic and completely integrable 
systems, the behaviors of power spectrum in the mixed 
regime between integrability and chaos is interesting. 
Interacting trapped bosons is a very complex system 
and due to the existence of two energy scales it nicely 
describes chaos to order transition with increase in 
number of energy levels. Our observation of Shnirelman 
peak strongly proves the earlier statement of Asaga 
[1T( . Our results nicely demonstrate how the degenerate 
single particle states of the pure harmonic trap are 
lifted gradually by increasing the effective interatomic 
interaction. Our findings are quite different from the 
results seen in atomic nuclei, atoms and molecules [l5j . 
Interacting trapped bosons is a very special and very 
complex system where the low-lying collective excita- 
tions are strongly influenced by interatomic interaction 
and shows level-repulsion. It is also spatially inhomo- 
geneous and the high-lying levels are of single particle 
nature and have regular features . For the dilute 
interacting Bose gas, it is also possible to calculate 
a large number of energy levels with high statistical 
precision. They can also be measured experimentally. 
The corresponding level fluctuation shows a transition 
from close to Wigner to Poisson with increase in energy 
levels showing it does not follow GOE predictions and 
we need a modified GOE which combines uniform, 
GOE and Poisson |47j. We observe the existence of 

The parameter 



i power law in the energy spectrum. 



a measures the fluctuation properties of the quantum 
system through the power spectrum. As the interacting 
trapped bosons arc interesting in the connection of 
recent experiments of BEC, our system is generic and it 
confirms that noise is ubiquitous in nature. However 
some open questions viz. how the spectral distribution 
will change with the attractive interactions to study 
the dynamical behavior of energy spectrum, still remain . 
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FIG. 2. (color online) Fig. l(a)-(c): Plot of fluctuation S n 
statistics for different energy levels. Fig. l(d)-(f): Plot of 
average power spectrum. 




FIG. 3. (color online) Histogram plot of P(s) distribution vs 
s for lowest 100 levels with different number of bosons in the 
trap. Blue lines are Wigner distribution. 
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FIG. 4. (color online) Level spacing distribution P(s) for the 
lowest 100 levels with 3 atoms in the trap. Fig. 4 (a)-(d) 
shows the Shnirelman peak in finer details. 




FIG. 5. (color online) Plot of integral level spacing distribu- 
tion I(s) againt ln(s) 
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FIG. 6. (color online) Histogram plot of P(s) distribution vs s 
for different the number of bosons (N) for energy-levels 3900 
to 4000. Green lines are Poisson distribution and Blue lines 
are Wigner distribution. 




FIG. 7. (color online) Spectral average < As(L) > computed 
for the Hamiltonian (1) with different number of interacting 
bosons (N) in the external trap vs L, (a) for lowest 100 energy 
levels and (b) for energy levels between 3900 and 4000. 



